/**
 * @file BDFs.h
 * @brief 
 * @author Linhuo (1066117119@qq.com)
 * @version 1.0
 * @date 2021-04-26
 * 
 * @copyright Copyright (c) 2021  Linhuo
 * 
 * @par 修改日志:
 * <table>
 * <tr><th>Date       <th>Version <th>Author  <th>Description
 * <tr><td>2021-04-26 <td>1.0     <td>Linhuo     <td>内容
 * </table>
 */

#ifndef _solution_BDFs_
#define _solution_BDFs_

#include "TimeIntegrator.h"

class BDF: public TimeIntegrator
{
private:
    Eigen::VectorXd alpha;
    Eigen::VectorXd beta;

public:
    void load_scheme(std::string _name);
    void One_step_forward();
};



void BDF::One_step_forward()
{   
    Eigen::MatrixXd I(6,6);
    I.setIdentity();
    Eigen::VectorXd u_new(6);
    Eigen::VectorXd _remain(6);
    Eigen::VectorXd u_temp_1(6);
    Eigen::VectorXd u_temp_2(6);
    u_new << 0,0,0,0,0,0;
    _remain << 0,0,0,0,0,0;
    u_temp_2 = u.col(s-1);
    for (int j = 0; j < s; j++)
	    _remain -= alpha(s-j-1) * u.col(j) ;
    Eigen::MatrixXd JF;
    JF.resize(6,6);
    Eigen::VectorXd b1;
    b1.resize(6);
    Eigen::VectorXd dif_u;
    dif_u.resize(6);
    dif_u << 100, 0, 0, 0, 0, 0;

    for(int i = 0; i < 5 && dif_u.norm() > 1e-5 ; i++)
    {   
        u_temp_1 = u_temp_2 ;
        JF = dt * beta(0) * Jfunc(u_temp_1) -I;
        b1 = u_temp_1 - dt * beta(0) * func(u_temp_1) - _remain;
        dif_u = (JF).partialPivLu().solve(b1); 
        u_temp_2 = dif_u + u_temp_1;
    }
    u_new = u_temp_2;

    for (int j = 0; j < s - 1; j ++)
    {
        u.col(j) = u.col(j+1);
        fu.col(j) = fu.col(j+1);
    };
    u.col(s-1) = u_new;
    fu.col(s-1) = func(u_new);
};


void BDF::load_scheme(std::string _name)
{
    bool is_found = false;
    std::ifstream scheme("lmm_scheme");
    while (!scheme.eof() && !is_found)
    {
	std::string tmp;
	std::getline(scheme, tmp);
	scheme >> name;
	if (name == _name)
	{
	    is_found = true;
	    scheme >> p;
        scheme >> s;
        scheme >> dt;
	    std::getline(scheme, tmp);
	    std::getline(scheme, tmp);
	    alpha.resize(s);
	    beta.resize(1);
        double alpha_Devisor=0,beta_Devisor=0;
        scheme >> alpha_Devisor;
	    for (int i = 0; i < s; i++)
	    {
		scheme >> alpha(i);
        alpha(i) /= (alpha_Devisor*1.0);
	    }
	    std::getline(scheme, tmp);
	    std::getline(scheme, tmp);
        scheme >> beta_Devisor;
		scheme >> beta(0);
        beta(0) /= (beta_Devisor*1.0);
	}
	else
	{
	    for (int i = 0; i < 6; i++)
		std::getline(scheme, tmp);
	}
    }
    if (!is_found)
    {
	std::cout << "scheme not found!" << std::endl;
	std::exit(-1);
    }
};


#else
//DO NOTHING
#endif